Ab initio variational solution for lattice thermal conductivity in diamond 



Giorgia Fugallo, Michele Lazzeri, Lorenzo Paulatto and Francesco Mauri 
IMPMC, Universite Pierre et Marie Curie, CNRS, 
4 place Jussieu, F -15252 Paris, France 



We present a first-principles theoretical approach for evaluating the lattice thermal conductivity 
based on the exact solution of the Boltzmann transport equation. We use the variational principle 
and the conjugate gradient scheme, which provide us with an algorithm faster than the one previously 
used in literature and able to always convergence to the exact solution. Three-phonon normal and 
umklapp collision, isotope scattering and border effects are rigorously treated in the calculation. 
Good agreement with experimental data for diamond is found. Moreover we show that that by 
growing more enriched diamond samples it is possible to achieve values of thermal conductivity up 
to three times larger than the commonly observed in isotopically enriched diamond samples with 
99.93% C 12 and 0.07 C 13 . 
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I. INTRODUCTION 

Thermal conductivity is one of the most important 
parameter used to characterize transport phenomena in 
solid state systems. A predictive theory for evaluating 
thermal conductivity is essential for the design of new 
materials for efficient thermoelectric refrigeration and 
power generation^ and it could help in understanding 
heat dissipation in micro- and nano-electronics devices!^ 

When heat is mostly carried by lattice vibrations, such 
as in semiconductors and insulators, a correct theoretical 
prediction of thermal transport properties cannot leave 
aside an accurate description of the phonon-phonon 
interactions and lifetimes. These quantities are related 
to second and third order derivatives of the ground 
state energy with respect to atomic displacements. 
Specifically the harmonic interatomic force constants 
determine phonon frequencies, group velocities and 
phonon populations while the anharmonic interatomic 
force constants determine phonon scattering rates and 
linewidths. 

A first microscopic description of the thermal con- 
ductivity in semiconductors and insulators has been 
formulated in 1929 by Peierls and it has became known 
as Boltzmann Transport Equation (BTE) . This equation 
involves the unknown perturbed population of a phonon 
mode and it describes how the perturbation due to a 
gradient of temperature is balanced by the change in the 
phonon population due to the scattering processes. A 
good predictive theory requires then a good knowledge 
of i) the harmonic and anharmonic inter-atomic force 
constants (IFCs) and ii) the perturbed phonon popula- 
tion obtained as solution of the BTE. 
Both these issues have non trivial solutions. The first 
issue can be addressed in the framework of Density 
Functional Perturbation Theory (DFPT) evaluating 
the interatomic force constants fully ab initio using the 
"2n+l" theorerrPHsl. An efficient implementation of 
this method extended for metallic systems exist in the 
Quantum EPRESSO package 6 for zone-centered mode^. 



A generalization for metallic systems and arbitrary 
phonons has recently been developed and implemented 
in the the Quantum-ESPRESSO packaged 
The second issue, lying in solving the BTE equation 
exactly, is due to the complexity of the scattering term. 
The change in the phonon population numbers of each 
single state involved in the scattering term depend in 
turn, on the change in the occupation number of the 
other states involved. 

Several theoretical studies instead of attempting to solve 
the BTE employ a common approximation, namely 
the single mode phonon relaxation time approximation 
(SMASHES This approximation describes rigorously 
the depopulation of the phonon states but not the 
corresponding repopulation, which is assumed to have 
no memory of the initial phonon distribution. The 
momentum conserving character of the normal (N) 
processes gives then rise to a conceptual inadequacy of 
the SMA description and its use becomes questionable 
in particular in the range of low temperatures where the 
umklapp (U) processes are frozen out and N processes 
lead the phonon relaxation^. 

Improved approxima te te chniques involve the use of a 
variational procedur e! 12 * 13 !. In such a kind of approach, 
originally introduced by Hamilton and ParrotlP^j, the 
thermal conductivity is found by variationally op- 
timizing a trial function adopted for describing the 
non-equilibrium phonon distribution function. Unfortu- 
nately the less the system is symmetric and isotropic the 
more the result and the accuracy will be affected by the 
form adopted for the trial function. 
A first approach to solve exactly the linearized BTE has 
been introduced in the 90s by Omini and Sparavigna^. 
The numerical solution evaluated on a discrete grid is 
obtained via a self-consistent iterative procedure, but 
as indicated by the authors^ there is not a general 
proof that convergence will always be obtained with 
this approach. Nevertheless until now the Omini Spar- 
avigna iterative procedure had represented the unique 
numerically exact method used to solve the BTE and 
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evaluating the thermal conductivity^"^. 

In this paper we present a new numerical approach for 
solving exactly the linearized BTE. This method joins 
together the variational principle and the resolution on a 
discrete grid. More specifically by using the variational 
principle and the conjugate gradient method we present 
a stable algorithm, faster than the one previously pro- 
posed and able to always converge to the exact solution. 
As a first application we use this algorithm for studying 
the lattice thermal conductivity in naturally occurring 
and isotopically enriched diamond. We find a good agree- 
ment with experiments and we are able to give for the 
first time a theoretical limit based on the exact solution 
of the BTE of the maximum lattice thermal conductivity 
reachable in isotopically pure diamond samples. 



II. BOLTZMANN TRASPORT EQUATION 

When a gradient of temperature VT is established in 
a system a subsequent heat flux will start propagating 
in the medium. Without loss of generality we assume 
that the gradient of temperature is along the direction x. 
The flux of heat, collinear to the temperature gradient, 
can be written in terms of phonon energies huiqj phonon 
group velocities c qj in the x direction and the perturbed 
phonon population riqj. 



qj °qj ' 'qj 



or 

" dx 



(1) 



in the l.h.s uj^j is the pulsation of the phonon mode with 
vector q and branch index j, Q is the volume of the unit 
cell and the sum runs over a uniform mesh of No q points. 
While in the r.h.s. k is the diagonal component of the 
thermal conductivity in the temperature-gradient direc- 
tion. Knowledge of the phonon perturbed population 
allows heat flux and subsequently thermal conductivity 
to be evaluated. 

Unlike phonon scattering by defects, impurities and 
boundaries the anharmonic scattering represents an in- 
trinsic resistive process and in high quality samples it 
dominates the behaviour of lattice thermal conductiv- 
ity at room temperature balancing the perturbation due 
to the gradient of temperature. The balance equation, 
namely the Boltzmann Transport Equation (BTE), for- 
mulated in 1929 by Peierls^l is: 



the equilibrium phonon population dn^/dT = dn^/dT 
where n^j = (e h ^^l Kl)T + while for the scattering 

term it can be expanded about its equilibrium value in 
terms of a first order perturbation: 
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The linearized BTE can be then written in the following 
fornPl 
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where the sum on q' and q" is performed in the Bril- 
louin Zone (BZ). The EX superscript denotes the exact 
solution on the BTE, to be distinguished to the approx- 
imated solutions we will discuss later. 
In Eq. [4] the anharmonic scattering processes as well as 
the scattering with the isotopic impurities and the bor- 
der effect are considered. More specifically (see Fig{T]) 

P^j q, ., is the scattering rate at the equilibrium of a pro- 
cesses where a phonon mode qj scatters by absorbing an- 
other mode q'j' to generate a third phonon mode q"j" 
. While P^j 3 q 3 is the scattering rate at the equilib- 
rium of a process where a phonon mode qj decays in two 
modes q'j' and q"j". 
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with the first term indicating the phonon diffusion due 
to the temperature gradient and the second term the 
scattering rate due to all the scattering processes. This 
equation has to be solved self consistently. In the gen- 
eral approach^, for small perturbation from the equilib- 
rium, the temperature gradient of the perturbed phonon 
population is replaced with the temperature gradient of 



Figure 1. Phonon scattering processes in a finite size anhar- 
monic crystal in presence of isotopic impurities. 

The two scattering rates have the forms: 
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with G the reciprocal lattice vectors. In order to evalu- 
ate them it is necessary to compute the third derivative 
of the total energy of the crystal £ tot ({w SQ ,(R;)}), 
with respect to the atomic displacement u SQ (R;), from 
the equilibrium position, of the s-th atom, along the a 
Cartesian coordinate in the crystal cell identified by the 
lattice vector R; : 



^ (3) (qj,qY,q"j") = 
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where £. cel1 is the energy per unit cell and the adimen- 
sional quantity X^j is defined by 
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z^j being the orthogonal phonon eigenmodes normalized 
on the unit cell and M s the atomic masses. 
The rate of the elastic scattering with isotopic impurities 
(see Fig|TJ) has the forrrPS 
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where gf = ^ Ms ^M^i^ is the average on the mass distri- 
bution of the atom of type s. In presence of two isotopes 
M s and M s > it can be written in terms of the concentra- 
tion e and mass change AAjf s = M s i — M s : 



<?! = e(l-e; 
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with (M s ) = M s + eAM s . 

Eventually, in system of finite size, P^° describes the 
absorption of phonon on the border (see Fig[I]): 



(11) 



where L is the Casimir length of the sample and F a 
correction factor depending on the width to length ratio 
of the sample. 

For the sake of clarity we will contract from now on the 
vector q and branch index j in one single mode index v. 
The BTE of Eq. [4] can be written as a linear system in a 
matrix form: 
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with the vector b v i = — c^Huj^i n u / (n v i +1) and the matrix 
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where we have used P^ — P^, v ,, from the detailed 
balance condition n u {n u i + l)(n„" + 1) = (n v + l)n v in v u 
(valid under the assumption hu = fiui' + fiu"). In this 
form the matrix is symmetric and positive semi-definite 
(see Appendix A for demonstrations) and it can be de- 
composed in A = A out + A m , where 
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tJ being the phonon relaxation time (see Appendix B). 
The A out diagonal matrix describes the depopulation of 
phonon states due to the scattering mechanisms while 
the A m matrix describes their repopulation due to the 
incoming output phonons. 

The solution of the linear system in Eq. [12] is obtained 
formally by inverting the matrix A. 
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and subsequently the thermal conductivity will be eval- 
uated as: 



k = Ab • f EX 
with A = h/N n k B T 2 . 



(17) 



III. SOLUTIONS OF THE BOLTZMANN 
TRANSPORT EQUATION 

The complexity of the BTE lies in the need of explicitly 
computing, storing and inverting the large matrix A. In 
the SMA the BTE is solved for the n v neglecting the role 
of the repopulation by means setting A ln to zero. 
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(18) 



Storing and inverting A out is trivial due to its diagonal 
form. The lattice thermal conductivity in SMA is then 



fc SMA = Ab . f ; 



SMA 
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Such approximation is exact if the repopulation loses 
memory of the initial phonon distribution and if it is 
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proportional to the equilibrium population of v. It re- 
mains anyway a good approximation if the repopulation 
is isotropic. An exact solution of Eq. [l2j that does not 
imply neither storing nor the explicit inversion of matrix 
A, has been proposed by Omini and Sparavigna (OS) 16 
by converging with respect to the iteration i the follow- 
ing: 
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with the iteration zero consisting in the SMA fo = f SMA . 
This procedure required, as for the SMA, only the trivial 
inversion of the diagonal matrix A out . Instead of stor- 
ing and inverting A, it just requires the evaluation of 
A ln fi, which is an operation computationally much less 
demanding, at each iteration i of the OS method. 
Once the convergence is obtained the thermal conductiv- 
ity is evaluated by : 



fc NV (f l ) = Ab-f l 



(21) 



From a mathematical point of views the OS iterative pro- 
cedure can be written as a geometric series: 
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that means that only if the absolute value of the ratio 
is smaller than one it is possible to reach a 
converged solution of the linear system in Eq. [12] . 
An alternative approach consists in observing that due 
to the properties of the matrix A (see Appendix A) it 
is possible to find the exact solution of the linearized 
BTE, Eq[l2j by using the variational principle. Indeed 
the solution of the BTE is t he ve ctor f EX which makes 
stationary the quadratic forrrpLUU 
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Since A is positive the stationary point is the global and 
single minimum of this functional. One can then define 
a variational conductivity functional: 



k v (f) = -2AJ(f) 



(24) 



that has the property fc v (f EX ) = k while any other value 
of fc v (f) underestimates k. In other words to find the 
minimum of the quadratic form it is equivalent to max- 
imize the thermal conductivity functional. As a conse- 
quence an error 8f = f — f EX results in an error in con- 
ductivity, linear in St if the functional is written in Eq. 
21 form, and quadratic if the variational form (Eq. 24) 
is used. 

In literature^, due to the complexity of the numeri- 
cal calculations the variational scheme was used together 
with a trial function for describing the non equilibrium 
phonon distribution function affecting then the accuracy 



of the final result with the form of the specific probe 
function chosen. In our scheme we avoid the use of trial 



function and we solve Eq. 12 on a grid (as in OS proce- 
dure) by using the conjugate gradient method^, as re- 
ported in Appendix C, to obtain the exact solution of the 
BTE equation. In order to speed up the convergence of 
the conjugate gradient we take advantage of the diagonal 
and dominant role of A out and we use a preconditioned 
conjugate gradient. Formally, this corresponds to use in 
the minimization the rescaled variable: 



f = VA° ut f 



(25) 



and then, with respect to this new variable, minimize the 
quadratic form ^(f) = ^(f) where: 
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Notice that f SMA = t>. The square root evaluation of 
A out is trivial due to its diagonal form. The computa- 
tional cost per iteration of the conjugate gradient scheme 
is equivalent to the OS one, but it always converges and 
requires a smaller number of iterations. 



IV. COMPUTATIONAL DETAILS 

In order to compute the thermal conductivity the only 
input required are the second and third order IFCs. 
Both of them were calculated by using the Quantum 
ESPRESSO package 6 following the method explained 
The first BZ is discretized into a uniform grid cen- 
tered at Gamma of q points. In such a way that if q and 
q' belong to the mesh also q ± q' belongs to the mesh, 
assuring a perfect momentum conservation. At any q 
the phonon frequencies are evaluated from the second 
order force constants and the phonon group velocities 
are computed from the derivative of the phonon disper- 
sion du/dq, using the Hellman-Feynmann theorem and 
obtaining the following velocity matrix directly from the 
Dynamical matrix D: 
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(29) 



In the non degenerate case c qj = Cjj while in the degen- 
erate one we use the phonon polarization vectors that 
diagonalize the matrix in the degenerate subspace. To 
compute the scattering rates, the BZ is again discretized 
into a grid of q' points centered in zero. The delta func- 
tion for the energy conservation is replaced by a Gaussian 
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It is important to notice that when the delta function 
is substituted with a Gaussian the detailed balance con- 
dition is only valid under approximation. This means 
that the OS definition of matrix A given irP^I and our 
definition, in Eq. 13 are not equivalent anymore. Our 



definition has the advantage to keep, for any finite a) 
in Eq. |30| the symmetric and non-negative character of 
the A matrix thanks to the symmetric definition of the 
scattering rate with the isotopic impurities given in Eq. 
[9] and the replacement of P„ ,v with P£, v „. 
For diamond calculations: a smearing a — 20 cm' 1 along 
the q' mesh of 30 x 30 x 30 has been found to lead to 
converged relaxation times and for the border effects: we 
used a Casim ir length L = 0.3 cm and a shape factor 

F = Q. j26T27l 

A norm conserving pseudopotentiaP^ with cutoff radii of 
1.2 a.u. and core correction has been used for C. The ex- 
change correlation energy is calculated in the framework 
of the Local Density Approximation (LDAp2). A plane- 
wave kinetic energy cutoff of 90 Ry and of 360 Ry for the 
charge density. We used a 8 x 8 x 8 Monkhorst-Pack of 
the BZ for the electronic k-point sampling. 
Anharmonic forces have been computed on a 4 x 4 x 4 
q-point phonon grid on the BZ, Fourier interpolated with 
a finer 30 x 30 x 30 mesh for the Boltzmann calculations. 
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Figure 2. (Color online). Lattice thermal conductivity of 
diamond at 100 K (top panel) and absolute error 5k (bottom 
panel) compared to the exact solution k for: the iterative 
scheme in the Omini Sparavigna standard definition (solid 
line), the iterative scheme in the variational definition given in 
Eq. [24] (dash-dotted line), and the conjugate gradient method 
with (dashed line) and without (dotted line) preconditioning, 
as a function of the order of iteration. 



V. RESULTS AND DISCUSSION 

In Fig[2]a comparison between the convergence trend 
obtained via the OS iteration scheme or the conjugate 
gradient is reported for the case of bulk diamond at 100 
K. The OS standard iterative scheme shows a numeri- 
cal instability after 77 iterations. This instability pre- 
vents the scheme to approach the exact solution k with 
a precision higher than ~ 300 W m' 1 K _1 . A higher 
precision is achievable with the Conj. Grad. already af- 
ter 4 iterations. As expected, if the variational definition 
of k (Eq. 24) is used in the OS iterative scheme, half 



the number of iterations are necessary to reach the same 
precision and the numerical instability appears after 41 
iterations. Moreover the convergence trend of the Conj. 
Grad. scheme without preconditioning is reported in the 
same graph to show how preconditioning is necessary to 
ensure a fast convergence. 

We have chosen for the comparison a temperature of 
100 K close to the maximum value of thermal conduc tiv- 
ity obtainable in finite size diamond samples ^ * 30 * 31 !. In 
this range of temperatures, where the U processes are a 
few and the border effects are not dominant, it is impor- 
tant to have a stable algorithm able to well characterize 
the few scattering processes that drive the lattice ther- 
mal conductivity in order to obtain the correct result. 

Fig. [3^i compares the lattice thermal conductivity of 
isotopically enriched and naturally occurring diamond 
obtained by solving exactly the BTE equation with the 
experimental results as a function of temperatures. The 
circles^ and squares^pointers represent the measured 
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Figure 3. (Color online). Lattice thermal conductivity 
of isotopically enriched and naturally occurring diamond 
as a function of temperatures. (Top) Experimental values 
(circled^, squares^ and triangles^) at different C 13 percent- 
ages (0.07%, 0.1% and 1.1%) are compared with the results 
of our ab initio calculations (solid, dash-dotted and dash with 
two dots lines). As indication of the theoretical limit a dashed 
line for the case in total absence of C 13 is reported. (Bottom) 
Ratio between the exact and the SMA solution as a function 
of temperature at different C 13 percentages. 



6 



values for isotopically enriched diamond with 99.93% 
C 12 , 0.07%C 13 and 99.9% C 12 , 0.1% C 13 respectively. 
While the triangles^ represent naturally occurring dia- 
mond with 98.9% C 12 and 1.1% C 13 . The results are 
in very good agreement. In the same picture is also in- 
dicated with a dashed line the thermal conductivity in 
total absence of C 13 . This value gives a theoretical limit 
of the maximum lattice thermal conductivity reachable. 
In the picture it is easy to notice that at 100 K, where 
the lattice thermal conductivity takes its maximum val- 
ues, & % cl3 — 3fco.o7% c i3 ■ This means that there is still 
a significant increment in thermal conductivity achiev- 
able by growing more enriched diamond samples. As the 
temperature increases, the values for the naturally occur- 
ring and isotopically pure samples become indistinguish- 
able. This is due to the U scattering becoming stronger 
and consequently driving the thermal conductivity as the 
temperature increases. For temperatures lower than 100 
K the border effects become dominant. 
FigjS]} shows the ratio between the thermal conductivity 
obtained by solving exactly BTE equation and by using 
the SMA as a function of temperatures (for T > 100 K). 
The lower the temperature and the less the C 13 abun- 
dance the bigger becomes the ratio between the exact 
solution and the SMA solution. In other words, the less 
are the events of scattering that do not conserve the mo- 
mentum (i.e. Umklapp, isotopes and border scattering) 
the less the SMA is able to give a good description of the 
process. In Figj3] is shown also the case with 99.995% 
C 12 and 0.005% C 13 as a further indication of how even 
small changes in the sample enrichment could give rise to 
sensible differences in the thermal transport properties of 
the material. 

In Fig|4] this last concept is more enlightened. Dia- 
mond thermal conductivity is represented as a function 
of isotopic presence for two different temperature 100 K 
and 300 K At T = 100 K the range of thermal conduc- 
tivity explores by changing the percentage of C 13 from 
to 1% spans one order of magnitude while at 300 K 
the differences between the different isotopic percentage 
becomes smaller with a ratio between the two extreme 
cases reduced of a factor 1.5. 



VI. CONCLUSIONS 

In this paper we have presented a new numerical ap- 
proach for solving exactly the linearized BTE. We have 
shown how joining the variational principle approach and 
the resolution on a grid it is possible to converge to the 
exact solution. Moreover the preconditioned conjugate 
gradient scheme with the line minimization assures a sig- 
nificant faster convergence than the method previously 
proposed by Omini and SparavigngP^ with an equivalent 
computational cost per iteration. 

As first application of our method we have then 
computed the lattice thermal conductivity of isotopically 
enriched and naturally occurring diamond by evaluating 




Figure 4. (Color online) Diamond lattice thermal conduc- 
tivity as a function of isotopic presence at 100 K (left) and 
300 K (right). The solid lines join the values obtained solv- 
ing the BTE while the dashed lines join the SMA solutions. 
The value of the thermal conductivity in absolute absence of 
C 13 in the two temperature cases is represented by a dotted 
orizontal line. 



the harmonic and anharmonic IFCs fully ab initio 
in the framework of DFPT using a recent general 
implementation of the "2n+l" theorem in the Quantum 
ESPRESSO package combined with an exact solution of 
the linearized phonon BTE. 

In agre emen t with what previously shown in 
literatur e! 20 * 26 ! we have demonstrated the severe in- 
adequacy of the commonly used SMA in the range 
of temperature [100 — 300] K for isotopically enriched 
diamond samples. In this range of temperatures, the 
lattice thermal conductivity shows a high sensitivity to 
the isotopic enrichment and our calculations suggest 
that by growing more enriched diamond samples it is 
possible to achieve values of thermal conductivity up 
to three times larger than the commonly observed in 
isotopically enriched diamond samples with 99.93% C 12 
and 0.07 C 13 . 
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Appendix A: Properties of matrix A 

It is easy to prove that the matrix A is symmetric 
A v y - A v t, v = by using the properties: P"y = P%, v 
and Pjffi = Pffi in the definition of A v y given in Eqjl3] 
It is also possible to prove that it is positive semi-definite. 
In order to show that, the matrix A in Eq |13| can be 
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written as: 

A = J2 P»,>»'y + E P vpV vy + J2 P» e V» (Al) 



is the relaxation time due to the anharmonic scattering 
processes with Tqj half width at half maximum of the 
corresponding phonon broadening; while 



where Dj^ v , is a matrix with all the element equal to zero 
apart those involving the triplets {v, v' , v"} 




whose eigenvalues are: 0, and 3. 

For the part representing the elastic scattering with the 
isotopes: T) v y is a matrix with all the elements equal to 
zero apart those involving the couples -{V, z/} 



D, 



v v 

1 -1 
-1 1 



(A3) 



with eigenvalues and 2. 

Finally for the border effect, D„ is a matrix with all the 
elements equal to zero apart those involving {v, u} 



D, 
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(A4) 



Since P" v i, P„ a °* and P„ c , are non negative then the 
total matrix is positive semi-definite because sum of pos- 
itive semi-definite matrices. 



Appendix B: Phonon relaxation times 
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is the relaxation time due to the border effects and 
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the relaxation time associated to the elastic scattering 
with isotopic impurities. 



Appendix C: Conjugate gradient method 



The conjugate gradient minimization 25 of Eq 23 or Eq. 
[26] requires the evaluation of the gradient g^ = Af; — b 
and a line minimization. Since the form is quadratic the 
line minimization can be done analytically and exactly. 
Moreover the information required by the line minimiza- 
tion at iteration i can be recycled to compute the gradient 
at the next iteration i + 1. Starting with an the initial 
vector f = f SMA , initial gradient g = Af - f SMA and 
letting ho = —go, the conjugate gradient method can be 
summarized with the recurrence: 



When different events of scattering are present such as 
anharmonic scattering, scattering with isotopic impurites 
and border effects the total phonon relaxation time 
is expressed by the Matthiessen's rule as: 



(fir 1 



(Bl) 



where: 



U = Ah; 

H+l — H ~~ l — n i 

h, ■ t, 



gi+i 
h i+ i 



Si+l 



hi ■ U 

Si+l ' Si+l ■ 



Si ' Si 



(CI) 
(C2) 

(C3) 
(C4) 



n ivo , ., ... 



2(nq/ 



(1 + fiq/ 



lq"j")S(Hu)qj — tkjjqiji 



(B2) 



where is the search direction and tj is an auxiliary 
vector. Notice that each iteration requires only one ap- 
plication of the matrix A on the vector h; as in the OS 
method. This is the computationally more demanding 
part of the conjugate gradient step. 
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